rm(list=ls(all=TRUE))
results <- read.csv("Data/Urchin_ClimateChange_Data.csv", header=T, na="NA")
library("lsmeans") #post-hoc tests
library("effects") #plot effects of modeling
library("lmtest") #linear mixed modeling
library("lme4") #linear mixed modeling
library("lmerTest") #calculate p-values in models
library("emmeans") #post-hoc tests
library("cowplot") #arrange plots
library("MuMIn")
library("car") #levenes tests, Anova, qqPlots,
library("tidyverse")
library("stats")
library("plyr") #splitting, applying, and combining data
library("dplyr")
library("onewaytests") #allows for Welch test with unequal variances (spine lengths)
library("stats")
library("effsize")
library("FSA")
library("rcompanion")
library("broom") #turns model output into data frames
library("plotrix") #Calculates standard error of the mean
##Initial Collection = Day -24
InitialDiameter <-
#seperate out the initial sizes of urchins on day -24, when desired conditions were met.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "-24") %>%
#sizes on day -24
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
InitialMean <- mean(InitialDiameter$Diameter)
InitialSE <- std.error(InitialDiameter$Diameter)
InitialMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=InitialDiameter)
anova(InitialMod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.0005035 0.0005035 1 19 0.1144 0.7389
## pH 0.0102269 0.0102269 1 19 2.3232 0.1439
## Temperature:pH 0.0003840 0.0003840 1 19 0.0872 0.7709
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(InitialMod))
## [1] 27 4
shapiro.test(residuals(InitialMod)) #no pass
##
## Shapiro-Wilk normality test
##
## data: residuals(InitialMod)
## W = 0.52108, p-value = 4.878e-11
hist(residuals(InitialMod)) #eh
# 2. Equal variances
leveneTest(residuals(InitialMod)~InitialDiameter$Treatment) #No pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.9592 0.04307 *
## 42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(InitialMod),resid(InitialMod,type="pearson"),col="blue")
Assemble Data
##After acclimating before ramp up = Day-13
AcclimatedDiameter <-
#seperate out the initial sizes of urchins on day -24, when desired conditions were met.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "-13") %>%
#sizes on day -24
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
AcclimatedMean <- mean(AcclimatedDiameter$Diameter)
AcclimatedSE <- std.error(AcclimatedDiameter$Diameter)
Results
AcclimatedMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=AcclimatedDiameter)
anova(AcclimatedMod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.015013 0.015013 1 19 0.2154 0.64782
## pH 0.254996 0.254996 1 19 3.6588 0.07097 .
## Temperature:pH 0.177404 0.177404 1 19 2.5455 0.12711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(AcclimatedMod)) #not super normal but... see below
## [1] 42 2
shapiro.test(residuals(AcclimatedMod)) #passes
##
## Shapiro-Wilk normality test
##
## data: residuals(AcclimatedMod)
## W = 0.9577, p-value = 0.09326
hist(residuals(AcclimatedMod)) #looks normalish
# 2. Equal variances
leveneTest(residuals(AcclimatedMod)~InitialDiameter$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.3961 0.2573
## 42
plot(fitted(AcclimatedMod),resid(AcclimatedMod,type="pearson"),col="blue")
Assemble Data
##After ramp up and desired conditions are reached to begin experiment = Day 1
Day1Diameter <-
#seperate out the initial sizes of urchins on day 1, when desired conditions were met.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "1") %>%
#sizes on day 1
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
Day1Mean <- mean(Day1Diameter$Diameter)
Day1SE <- std.error(Day1Diameter$Diameter)
Results
Day1Mod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=Day1Diameter)
anova(Day1Mod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.04177 0.04177 1 19 0.1922 0.66606
## pH 0.38094 0.38094 1 19 1.7525 0.20128
## Temperature:pH 0.89969 0.89969 1 19 4.1389 0.05611 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(AcclimatedMod)) #not super normal but... see below
## [1] 42 2
shapiro.test(residuals(AcclimatedMod)) #passes
##
## Shapiro-Wilk normality test
##
## data: residuals(AcclimatedMod)
## W = 0.9577, p-value = 0.09326
hist(residuals(AcclimatedMod)) #looks normalish
# 2. Equal variances
leveneTest(residuals(AcclimatedMod)~InitialDiameter$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.3961 0.2573
## 42
plot(fitted(AcclimatedMod),resid(AcclimatedMod,type="pearson"),col="blue")
Subset data frame to obtain growth measurements.
### GROWTH MODEL 1: using day 1 as the first official day of experiment when desired environmental conditions were met (24 days after getting urchins. In that 24 days, they were acclimated and then conditions ramped up)
Initial <-
#seperate out the initial sizes of urchins on day 1, when desired conditions were met to standardize growth to this measurement.
results %>%
select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "1") %>%
#sizes on day 1
select("Temperature", "pH", "Diameter", "TankID")
mean(Initial$Diameter)
## [1] 16.03
std.error(Initial$Diameter)
## [1] 0.2550688
Final <-
#seperate out the final sizes of urchins on day 126
results %>%
select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column
filter(Day == "126") %>%
#filter for sizes on day 126
select("Temperature", "pH", "Diameter","TankID")
Growth <-
#create column that is growth
bind_cols(Initial, Final) %>%
#binds initial (Diameter) and final (Diameter1) sizes to calculate growth
mutate(Growth = ((Diameter1 - Diameter)/Diameter*100)) %>%
#add column for growth calculation
select("Growth")
resultsGrow <-
#table of initial and final size data
results %>%
select("Day", "Treatment", "Temperature", "pH", "TankID", "Diam1", "Diam2") %>%
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
filter(Day == "126")
resultsGrow <- bind_cols(resultsGrow, Growth)
#combines table of initial and final with growth information
#Figure out means of each treatment group
resultsGrow %>%
dplyr::group_by(Temperature, pH) %>%
dplyr::summarise(mean = mean(Growth), s.e. = se(Growth))
## # A tibble: 4 x 4
## # Groups: Temperature [2]
## Temperature pH mean s.e.
## <fct> <fct> <dbl> <dbl>
## 1 ambient acidified 327. 16.3
## 2 ambient ambient 323. 11.0
## 3 heated acidified 352. 15.0
## 4 heated ambient 376. 7.15
GrowthMean <- mean(resultsGrow$Diameter)
#total grpwth mean pooled for all treatments - not that necessary
Day1SE <- std.error(Day1Diameter$Diameter)
#same as above for SE
Results
#LMM for growth:
modGrow <-
resultsGrow %>%
lmer(Growth~ Temperature * pH + (1|TankID), data=.)
#Need to decide on model type that we want to use.
summary(modGrow)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 417.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.90063 -0.32406 0.08141 0.32903 1.94882
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 1929.2 43.92
## Residual 283.1 16.83
## Number of obs: 46, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 326.595 18.578 19.000 17.580 3.28e-13
## Temperatureheated 25.439 26.273 19.000 0.968 0.345
## pHambient -3.424 26.273 19.000 -0.130 0.898
## Temperatureheated:pHambient 27.731 38.073 19.000 0.728 0.475
##
## (Intercept) ***
## Temperatureheated
## pHambient
## Temperatureheated:pHambient
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707
## pHambient -0.707 0.500
## Tmprtrhtd:H 0.488 -0.690 -0.690
anova(modGrow, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 1169.41 1169.41 1 19 4.1304 0.05634 .
## pH 74.92 74.92 1 19 0.2646 0.61290
## Temperature:pH 150.20 150.20 1 19 0.5305 0.47527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of normality of residuals and homogeneity of variance. Assumptions are met.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow))
## [1] 43 20
shapiro.test(residuals(modGrow)) #passes
##
## Shapiro-Wilk normality test
##
## data: residuals(modGrow)
## W = 0.95491, p-value = 0.07263
hist(residuals(modGrow)) #looks normal
# 2. Equal variances
leveneTest(residuals(modGrow)~resultsGrow$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.1906 0.1033
## 42
plot(fitted(modGrow),resid(modGrow,type="pearson"),col="blue")
Display exploratory interaction plot.
#Interaction plot shows slight interaction
interaction.plot(resultsGrow$Temperature,resultsGrow$pH,resultsGrow$Growth,
type = c("l","p", "b", "o","c"), legend = TRUE,
trace.label = deparse(substitute(pH)),
fixed = FALSE,
xlab = deparse(substitute(Temperature)),
ylab = deparse(substitute(Growth)),)
Model growth using the start as the day conditions began to ramp up.
### GROWTH MODEL 2: If run the same model using day -13 as the first technical day of the experiment. This is the day conditions began to ramp up, so urchins were experiencing gradual increase of stressors.
Initial1 <-
#seperate out the initial sizes of urchins on day -14, when conditions began to ramp up.
results %>%
select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column
filter(Day == "-13") %>%
#filter for day at -13 to use this as the intial size.
select("Temperature", "pH", "Diameter", "TankID")
mean(Initial1$Diameter)
std.error(Initial1$Diameter)
Final1 <-
#seperate out the final sizes of urchins on day 126
results %>%
select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column
filter(Day == "126") %>%
#filter for size on last day of experiment
select("Temperature", "pH", "Diameter","TankID")
Growth1 <-
#create column that is growth
bind_cols(Initial1, Final1) %>%
#binds initial (Diameter) and final (Diameter1) sizes to calculate growth
mutate(Growth = ((Diameter1 - Diameter)/Diameter*100)) %>%
#add column for growth calculation
select("Growth")
resultsGrow1 <-
#table of initial and final size data
results %>%
select("Day", "Treatment", "Temperature", "pH", "TankID", "Diam1", "Diam2") %>%
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
filter(Day == "126")
resultsGrow1 <- bind_cols(resultsGrow1, Growth1)
#combines table of initial and final with growth information
#Figure out means of each treatment group
resultsGrow1 %>%
dplyr::group_by(Temperature, pH) %>%
dplyr::summarise(mean = mean(Growth), s.e. = se(Growth))
Analyze with linear mixed model.
#LMM for growth:
modGrow1 <-
resultsGrow1 %>%
lmer(Growth~ Temperature * pH + (1|TankID), data=.)
summary(modGrow1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 451.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.97009 -0.36352 -0.06635 0.38074 2.05538
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 4120 64.19
## Residual 667 25.83
## Number of obs: 46, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 556.515 27.244 19.000 20.427 2.17e-14
## Temperatureheated 56.177 38.529 19.000 1.458 0.161
## pHambient -3.451 38.529 19.000 -0.090 0.930
## Temperatureheated:pHambient 9.115 55.833 19.000 0.163 0.872
##
## (Intercept) ***
## Temperatureheated
## pHambient
## Temperatureheated:pHambient
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707
## pHambient -0.707 0.500
## Tmprtrhtd:H 0.488 -0.690 -0.690
anova(modGrow1, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 3141.63 3141.63 1 19 4.7100 0.04287 *
## pH 0.68 0.68 1 19 0.0010 0.97489
## Temperature:pH 17.78 17.78 1 19 0.0267 0.87205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow1))
## [1] 2 27
shapiro.test(residuals(modGrow1)) #passes
##
## Shapiro-Wilk normality test
##
## data: residuals(modGrow1)
## W = 0.95274, p-value = 0.05986
hist(residuals(modGrow1)) #looks normal
# 2. Equal variances
leveneTest(residuals(modGrow1)~resultsGrow1$Treatment) #doesn't pass...
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.9964 0.0413 *
## 42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modGrow1),resid(modGrow1,type="pearson"),col="blue")
Display exploratory interaction plot.
#Interaction plot shows interaction
interaction.plot(resultsGrow1$Temperature,resultsGrow1$pH,resultsGrow1$Growth,
type = c("l","p", "b", "o","c"), legend = TRUE,
trace.label = deparse(substitute(pH)),
fixed = FALSE,
xlab = deparse(substitute(Temperature)),
ylab = deparse(substitute(Growth)),)
Analysis of calcification at the tip of the spine:
Assemble data.
### Model 1: calcification at the tips of spines
resultsTip <-
#create data using only the SEM images of the tips of spines.
results %>%
select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>%
#selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
drop_na(RatioSEM) %>%
filter(Chosen == "yes", PartOfSpine == "Tip")
#filters only for tips and chosen side. No dusty images to be analyzed.
Analyse with linear mixed model:
#LMM for calcification at the tips of spines
modTip <-
resultsTip %>%
lmer(RatioSEM ~ Temperature * pH + (1|TankID), data = .)
## boundary (singular) fit: see ?isSingular
summary(modTip)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 64.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76413 -0.74507 -0.02052 0.63045 2.07630
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 0.0000 0.0000
## Residual 0.1356 0.3682
## Number of obs: 67, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.69118 0.08930 63.00000 18.938 <2e-16
## Temperatureheated -0.06662 0.12453 63.00000 -0.535 0.5945
## pHambient -0.07462 0.12453 63.00000 -0.599 0.5512
## Temperatureheated:pHambient 0.30557 0.18089 63.00000 1.689 0.0961
##
## (Intercept) ***
## Temperatureheated
## pHambient
## Temperatureheated:pHambient .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.717
## pHambient -0.717 0.514
## Tmprtrhtd:H 0.494 -0.688 -0.688
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(modTip, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.10158 0.10158 1 63 0.7492 0.39000
## pH 0.08185 0.08185 1 63 0.6038 0.44006
## Temperature:pH 0.38685 0.38685 1 63 2.8534 0.09612 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modTip)) #Normal
## [1] 38 8
shapiro.test(residuals(modTip)) #Pass
##
## Shapiro-Wilk normality test
##
## data: residuals(modTip)
## W = 0.97737, p-value = 0.2605
hist(residuals(modTip))
# 2. Equal variances
leveneTest(residuals(modTip)~resultsTip$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.4545 0.715
## 63
plot(fitted(modTip),resid(modTip,type="pearson"),col="blue")
Conduct post hoc analysis.
## POST HOC ANALYSIS
emmTip <- emmeans(modTip, ~Temperature*pH, adjust = "tukey") #Estimated marginal means (Least-squares means)
multcomp::cld(emmTip) #create a compact letter display of all pair-wise comparisons
## Temperature pH emmean SE df lower.CL upper.CL .group
## ambient ambient 1.62 0.0868 18.2 1.43 1.80 1
## heated acidified 1.62 0.0868 18.2 1.44 1.81 1
## ambient acidified 1.69 0.0900 18.2 1.50 1.88 1
## heated ambient 1.86 0.0993 16.6 1.65 2.07 1
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
pwpp(emmTip) #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.
Display interaction plot for exploration.
#Interaction plot for cacification at the tip of the spines - interaction
interaction.plot(resultsTip$Temperature,resultsTip$pH,resultsTip$RatioSEM,
type = c("l","p", "b", "o","c"), legend = TRUE,
trace.label = deparse(substitute(pH)),
fixed = FALSE,
xlab = deparse(substitute(Temperature)),
ylab = deparse(substitute(Growth)),)
Assemble data.
### Model 2: calcification in the base of the spines
resultsBase <-
#create data using only the SEM images of the base of spines.
results %>%
select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>%
#selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
drop_na(RatioSEM) %>%
filter(Chosen == "yes", PartOfSpine == "Base")
#filters only for tips and chosen side. No dusty images to be analyzed.
Results
#LMM for calcification at the tips of spines
modBase <-
resultsBase %>%
lmer(RatioSEM~ Temperature * pH + (1|TankID), data=.)
summary(modBase)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 98.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5535 -0.6315 -0.1529 0.6918 1.9937
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 0.1303 0.3610
## Residual 0.1597 0.3996
## Number of obs: 68, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.1440 0.1749 18.8077 12.257 2.06e-10
## Temperatureheated 0.3566 0.2487 19.1707 1.434 0.1677
## pHambient 0.6524 0.2474 18.8077 2.637 0.0163
## Temperatureheated:pHambient -0.2224 0.3594 18.9804 -0.619 0.5434
##
## (Intercept) ***
## Temperatureheated
## pHambient *
## Temperatureheated:pHambient
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.703
## pHambient -0.707 0.497
## Tmprtrhtd:H 0.487 -0.692 -0.688
anova(modBase, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.30743 0.30743 1 18.992 1.9251 0.181362
## pH 1.48873 1.48873 1 18.960 9.3223 0.006553 **
## Temperature:pH 0.06114 0.06114 1 18.980 0.3829 0.543424
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modBase)) #Normal
## [1] 33 37
shapiro.test(residuals(modBase)) # Pass
##
## Shapiro-Wilk normality test
##
## data: residuals(modBase)
## W = 0.97304, p-value = 0.1472
hist(residuals(modBase))
# 2. Equal variances
leveneTest(residuals(modBase)~resultsBase$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.068 0.369
## 64
plot(fitted(modBase),resid(modBase,type="pearson"),col="blue")
Post hoc analysis
## POST HOC ANALYSIS
emmBase <- emmeans(modBase, ~Temperature*pH, adjust = "tukey")
#Estimated marginal means (Least-squares means)
multcomp::cld(emmBase)
## Temperature pH emmean SE df lower.CL upper.CL .group
## ambient acidified 2.14 0.175 18.8 1.78 2.51 1
## heated acidified 2.50 0.177 19.5 2.13 2.87 12
## ambient ambient 2.80 0.175 18.8 2.43 3.16 12
## heated ambient 2.93 0.192 18.8 2.53 3.33 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
#create a compact letter display of all pair-wise comparisons
pwpp(emmBase)
#Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.
Display interaction plot.
#Interaction plot for calcification at the base of the spines - no interaction
interaction.plot(resultsBase$Temperature,resultsBase$pH,resultsBase$RatioSEM,
type = c("l","p", "b", "o","c"), legend = TRUE,
trace.label = deparse(substitute(pH)),
fixed = FALSE,
xlab = deparse(substitute(Temperature)),
ylab = deparse(substitute(RatioSEM)),)
Assemble data.
## Number of loose spines counted at the bottom of tanks at the end of the experiment. These were either broken or shed. Note: these spines were not collected, just counted, so can not determine whether they were broken or shed.
resultsDropped <-
#create data using only the needed variables.
results %>%
select ("Treatment", "Temperature", "pH", "SpineCount", "TankID") %>%
drop_na(SpineCount)
Analyse with linear mixed effect model.
##LM of dropped spines - can't make it LMM because error: number of levels of each grouping factor must be < number of observations - only have one count for each TankID.
modDropped <-
resultsDropped %>%
lm(SpineCount~ Temperature * pH, data=.)
summary(modDropped)
##
## Call:
## lm(formula = SpineCount ~ Temperature * pH, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0000 -6.0833 -0.1667 0.4000 20.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.167 3.640 5.814 1.33e-05 ***
## Temperatureheated -1.167 5.148 -0.227 0.82315
## pHambient -21.000 5.148 -4.079 0.00064 ***
## Temperatureheated:pHambient 1.600 7.461 0.214 0.83248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.917 on 19 degrees of freedom
## Multiple R-squared: 0.6088, Adjusted R-squared: 0.547
## F-statistic: 9.855 on 3 and 19 DF, p-value: 0.0003899
anova(modDropped) #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Analysis of Variance Table
##
## Response: SpineCount
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 1 1.52 1.52 0.0192 0.8914
## pH 1 2345.78 2345.78 29.4995 3.062e-05 ***
## Temperature:pH 1 3.66 3.66 0.0460 0.8325
## Residuals 19 1510.87 79.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##ASSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modDropped)) #not normal
## 1458 1460
## 2 4
shapiro.test(residuals(modDropped)) #faaaiiilll
##
## Shapiro-Wilk normality test
##
## data: residuals(modDropped)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped))
# 2. Equal variances
leveneTest(residuals(modDropped)~resultsDropped$Treatment) # barely pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.8487 0.06483 .
## 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modDropped),resid(modDropped,type="pearson"),col="blue")
Analysis does not meet assumptions. Attempt a square root transformation.
##DOES NOT MEET ASSUMPTIONS so ATTEMPT AT Transformation:
resultsDropped$tdata<-(resultsDropped$SpineCount)^1/2
# transformation
modDropped1 <- lm(tdata~ Temperature * pH, data=resultsDropped)
# run lm model of transformed data
###ASSUMPTIONS for Transformation
# 1. Normality of residuals
qqPlot(residuals(modDropped1)) #not normal
## 1458 1460
## 2 4
shapiro.test(residuals(modDropped1)) #fail
##
## Shapiro-Wilk normality test
##
## data: residuals(modDropped1)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped1))
# 2. equal variances
bartlett.test(residuals(modDropped1)~resultsDropped$Treatment) #fail
##
## Bartlett test of homogeneity of variances
##
## data: residuals(modDropped1) by resultsDropped$Treatment
## Bartlett's K-squared = 43.282, df = 3, p-value = 2.144e-09
plot(fitted(modDropped1),resid(modDropped1,type="pearson"),col="blue")
summary(modDropped1)
##
## Call:
## lm(formula = tdata ~ Temperature * pH, data = resultsDropped)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0000 -3.0417 -0.0833 0.2000 10.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.5833 1.8202 5.814 1.33e-05 ***
## Temperatureheated -0.5833 2.5742 -0.227 0.82315
## pHambient -10.5000 2.5742 -4.079 0.00064 ***
## Temperatureheated:pHambient 0.8000 3.7304 0.214 0.83248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.459 on 19 degrees of freedom
## Multiple R-squared: 0.6088, Adjusted R-squared: 0.547
## F-statistic: 9.855 on 3 and 19 DF, p-value: 0.0003899
Anova(modDropped1, type=2) #test for significance of model
## Anova Table (Type II tests)
##
## Response: tdata
## Sum Sq Df F value Pr(>F)
## Temperature 0.23 1 0.0118 0.9146
## pH 586.44 1 29.4995 3.062e-05 ***
## Temperature:pH 0.91 1 0.0460 0.8325
## Residuals 377.72 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Transfomraiton was not successful, so we use a non parametric Kruskal Wallis Test.
### Transformations not work, so use non parametric:
kruskal.test(resultsDropped$SpineCount~resultsDropped$Treatment)# SIG*** (p=0.0005563)
##
## Kruskal-Wallis rank sum test
##
## data: resultsDropped$SpineCount by resultsDropped$Treatment
## Kruskal-Wallis chi-squared = 17.505, df = 3, p-value = 0.0005563
Conduct Dunn Post Hoc Test and display interaction plot.
### POST HOC Pairwise ######
dunn <-
resultsDropped %>%
dunnTest(SpineCount ~ Treatment,
method = "holm", kw = TRUE, data = .)
dunnph <- dunn$res
cldList(P.adj ~ Comparison, data = dunnph, threshold = 0.05)
## Group Letter MonoLetter
## 1 AMB a a
## 2 OA b b
## 3 T ac a c
## 4 T/OA bc bc
#Interaction plot for dropped spines - no interaction
interaction.plot(resultsDropped$Temperature,resultsDropped$pH,resultsDropped$SpineCount,
type = c("l","p", "b", "o","c"), legend = TRUE,
trace.label = deparse(substitute(pH)),
fixed = FALSE,
xlab = deparse(substitute(Temperature)),
ylab = deparse(substitute(SpineCount)),)